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Abstract 

Numerical solutions of a non-Fickian diffusion equation belonging to a hyperbolic 
type are presented in one space dimension. The Brownian particle modelled by this 
diffusion equation is subjected to a symmetric periodic potential whose spatial shape 
can be varied by a single parameter. We consider a numerical method which consists 
of applying Laplace transform in time; we then obtain an elliptic diffusion equation 
which is discretized using a finite difference method. We analyze some aspects of the 
convergence of the method. Numerical results for particle density, flux and mean- 
square-displacement (covering both inertial and diffusive regimes) are presented. 

keywords: numerical methods, Laplace transform, telegraph equation, periodic 
potential, non-Fickian diffusion 

1 Introduction 

In this paper we shall present numerical solutions of a non-Fickian diffusion equation in 
the presence of a symmetric periodic potential in one space dimension. Let us briefly recall 
that the Fickian diffusion equation in the presence of a potential V(£) reads 



dn ,, x r^d 2 n.^ . 1 d 



where 7 is a friction parameter and D = ksT/m^ is the diffusion coefficient, m being the 
mass of the Brownian particle whose overdamped (diffusive) dynamics is well described 
by ([T]), ks is the Boltzmann's constant and T the temperature of the fluid. 
The equation of our study is 

l<9 2 n_ . On.. . ^d 2 n.^ . 1 d 
:(M + —(£,t) = D—(£,t) + 



(2) 



7 <9t 2 ' dr ' d£ 2 ' rwy d(, 

Both equations, (pQ) and ([2]) can be derived from an underlying kinetic equation e.g. the 
phase-space Kramers equation [S] 

dr m <9£ d£ dp ^ dp ^ 171 b 1 q^ 2 ' 
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where t) is the probability density function for the position component £ and mo- 

mentum component p of the Brownian particle. 

Equation ([2j) in the absence of a potential field is sometimes referred to the telegrapher 
equation although we shall call it a non-Fickian diffusion equation. We refer to [8] f° r a 
derivation of (|2]) from ([3]). It may be noted that for times larger than I/7 , the first term 
on the left hand side of ([2]) can be neglected and the Fickian regime is regained. Equation 
([T]) in the absence of a potential field leads to the well known result for the mean square 
displacement [IDJ 

< f( T ) >= 2Dt. (4) 

In the presence of a flexible symmetric potential, it was shown in [8] that < £ 2 (r) > does 
not necessarily behave linearly with time. Equation retains some short time inertial 
behaviour of a Brownian particle and at long time results in a diffusive behaviour. The 
velocity v = d^/dr of a Brownian particle is not well defined in the diffusive regime for 
which ([I]) is applicable. Since (J2]) is applicable in an inertial regime, the velocity can be 
calculated with Quite recently the instantaneous velocity of a Brownian particle has 
been experimentally investigated [11} [T2] [17] . This provides an additional motivation for 
studying ([2]). There is also a recent paper [1] which models transport of ions in insulating 
media through a non-Fickian diffusion equation of the type discussed in our work. In [3] 
the non-Fickian diffusion equation is referred to as a hyperbolic diffusion equation . 

To solve our problem we consider a numerical method based on space discretization 
and time Laplace transform. The latter is suitable for long times and also for solutions 
that are not necessarily smooth in time. It may be noted that iterative methods in time, 
including implicit methods such as the Crank- Nicolson [7J , which allows a choice of large 
time steps, usually take too long to compute the solution. 

The paper is organized as follows. In section 2 we present the model problem in 
dimensionless variables. In section 3 we describe a numerical method based on the time 
Laplace transform which is suitable for long time integration and also for solutions which 
are not very smooth. In section 4 the convergence properties of the algorithm are studied. 
In section 5 we present the behaviour of the solution to the non-Fickian diffusion equation, 
the flux and the mean square displacement. We conclude the paper, in section 6, with a 
summary and outlook. 



2 The model and physical quantities 

In our studies we consider three quantities of physical interest. These are the particle 
density n(£,r), the current density (flux) j(£,t) and the mean square displacement 
< £ 2 (t) >. The current density is not normally studied. However, since we are dealing 
with a non-Fickian diffusion equation we have decided to consider j(£,r) as well. For 
the Fickian case and in the absence of any potential, j(£,r) is related to n(£,r) through 
j = —D(dn/d£). This is not so in the non-Fickian case for which the relation between 
j(£,t) and n(£,r) is more involved. 

Let us consider the non-Fickian diffusion equations for particle density and the flux 

1 d 2 n dn 1 d .„ . d 2 n . . 

7 or 1 or rri"f a£ a£ 2 

1 dj „dn 1 „ 
J + -7T = -D^z Pn, 6 

with n(£, r) as the density of the Brownian particles. P is the force acting on the particle 
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due to the potential field V, i.e. 



dV_ 



We consider a symmetric periodic potential field, as previously studied in [5], [8] and 
[T3] . It reads 

F(£;a) = — ?— e^-1, (7) 
J (ia) 

where Jo is the Bessel function of the first kind and zero order and i is the imaginary 
number. In order to illustrate the flexible form of this single-parameter potential we have 
plotted, Figure dj the potential (J7|), for two values of the parameter, a = 1 and a = 16. 
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Figure 1: Potential field V(x;a). Left figure: a = 1; Right figure: a = 16. 

Our model consists of equations (J5]) and ([6]), and the potential field V(£;a) given by 
(|7|). For later purpose we introduce the following dimensionless parameters 



n 



n 



n 



i = T7, 



(8) 



where no is a reference particle density (concentration). The dimensionless forms of ([2j) 
and (El) can be written as 



d 2 n dn 

W + ~dt 

dj 



J + 



with 



dt 

P(x) 



d d 2 n 
dx dx 2 ' 



dV 



m^D-f 3 dx 



(9) 
(10) 

(11) 



3 Numerical method 

We consider equations ([9]) and (|10p with the following initial conditions 



n (x, 0) 
j(x,0) 



1 



1 

I7^ f 



_ X 2 /L 2 dn 



2 /L 2 



dt 



(x,0) = 0, 



L 2 



2x + jP(x) 



(12) 
(13) 
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where 

P{x) = 1=^, V(x;a) = ^_ e — - 1. 

The boundary conditions are given by 

lim n(x,t) = 0, lim n(x,t) = (14) 

x— >oo x— y— oo 

and 

lim j(x,t) = 0, lim j(x,t) = 0. (15) 

:r— >oo :r— oo 

In this section we describe a numerical method to solve the problem (f9|)- (|10p . Our 
approach can be separated in three steps. First, we apply the Laplace transform to ([9])- 
(|10p in order to remove the time dependent terms and we obtain an ordinary differential 
equation in x that also depends on the Laplace transform parameter s. Secondly, we solve 
the ordinary differential equation obtained using a finite difference scheme. Lastly, using a 
numerical inverse Laplace transform algorithm we obtain the final approximate solution. 

3.1 Spatial discretization 

Our numerical method is facilitated if we apply time Laplace transform to equation (|9|) 
and obtain the ordinary differential equation 

p±-\ s n-± {P n) = -(1 + s)n(x, 0), (16) 
aar ax 

where A s = s 2 + s, s is a complex variable and n is the Laplace transform of n defined by 

oo 



n(x,s) = / e s n(x,t)dt. 
Jo 

Now, assume we have a space discretization Xi = a + iAx, i = 0, . . . , N. Let rji(s), i = 
0, . . . , N represent the approximation of n (xi, s) in the Laplace transform domain. The 
outflow boundary is such that t]n( s ) = 0, for all s and N sufficiently large, which is 
according to the physical boundary condition. 

To derive the numerical method we consider central differences to approximate the 
first derivative and the second derivative of equation (TT6|) . We obtain, for a fixed s, the 
finite difference scheme given by 

^ l(S) = ^±^^ .x 9 9h{s)- ^^^=^=^ = -(l + sM Xi ,0), (17) 

for i = 1, . . . , N — 1, where Pi = P(xi). 
Therefore, we obtain the linear system 

K(s)r}(s) = b(s), (18) 

where K(s) = [Kij(s)] is a band matrix of size N — 1 x iV — 1 with bandwidth three and 
rj(s) = [rji (s) , . . . , t/at-i (s)] T . The matrix K(s) has entries of the form 

is ( \ 1 i 



Ax 2 2Ax ' 
2 
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and b (s) contains boundary conditions, being represented by 



6 (a) 



-(l + s)n(a?i,0) 
-(1 + s)n(x 2 ,0) 

-(1 + s)n(xjv_2,0) 
-(1 + s)n(xAr_i,0) 



+ 



-K li0 (s)r} (s) 






(20) 



-K N - liN (s)i] N (s) 

To compute the flux, we apply the Laplace transform to equation (fTUj) . that is, 
(1 + s)3 = -"^V^-^ + lP{x)n + j(x, 0), 



(21) 



where j is the Laplace transform of the flux j. The last step is to determine an approximate 
solution rji (t) and ji(t) of n(xi,t) and j(xi,t) respectively, which is obtained from rji (s) 
and ji (s) by using a Laplace inversion numerical method. 

3.2 Laplace transform inversion 

In this section, we determine an approximate solution r]i (t) from rji (s) by using a Laplace 
inversion numerical method. For the sake of clarity we omit the index i, denoting rji (s) 
by rj(s). 

A formally exact inverse Laplace transform of rj (s) into rj (t) is given through the 
Bromwich integral [14] 

-1 f-/3+ioo 

rj(t) = — e st ?j(s)ds, (22) 

where f3 is such that the contour of integration is to the right-hand side of any singularity 
of rj (s). However, for a numerical evaluation the above integral is first transformed to an 
equivalent form 



rj(t) 



1 



7T 



Re{^(s)e iw *} du, 



(23) 



where s = f3 + ico [H [TH [15]. The integral is now evaluated through the trapezoidal rule 
[TJE], with step size tt/T, and we obtain 



k=l 



_ . ifcvT \ ikTvt 

v ( P + ~y I e T 



Et, 



(24) 



for < t < 2T and where Et is the discretization error. It is known that the infinite 
series in this equation converges very slowly. To accelerate the convergence, we apply the 
quotient-difference algorithm, proposed in [2], and also used in [15], to calculate the series 
in (|24j) by the rational approximation in the form of a continued fraction. Under some 
conditions we can always associate a continued fraction to a given power series. 
We denote v (z) the continued fraction 



v (z) = d / (1 + d lZ / (1 + d 2 z/ (1 + • • • ))) 
associated to the power series in (|24p . For z = e 17Tt / T , 



v [z 



k=l 



iktr 



(25) 



(26) 
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and the coefficients di's of (|25p are obtained by recurrence relations from the coefficients 
Let the M-th partial fraction be denoted by v(z, M). Therefore 



v(z)=v(z,M)+EP, 
where Ep is the truncation error. Then 

rj (t) = ±e#Re {v (z, M) + Ep 1 } - E T . 

The approximation for fj (t) is denoted by n(t) and given by 

V (t) = ^ t Re{v(z,M)}. 

4 Convergence of the numerical method 

In this section we discuss the convergence of the numerical method chosen to compute an 
approximate solution to equation ([9]). Let us denote by E$ the error associated to the 
spatial discretization, that is, 

n(xi,s) = rji(s) + E s (xi,s). (27) 

The next errors come from the numerical inversion of Laplace transform, where the Laplace 
inverse transform of rji(s) is, as described in the previous section, the solution 

Vi(t) = ie^'Re^^MO + ^Cxi.t)} - E T ( Xi ,t), (28) 

where Et is the error associated with the trapezoidal approximation and Ep is the trun- 
cation error associated to the continued fraction. Note that for each Xi the algorithm 
chooses a Mj and therefore for each Xi we have a different value of the approximation of 
the continued fraction, v (z,Mi). Therefore from (|27p - (|28|) we have 



n( Xi ,t) = ^ t Re{v(z,M i ) + E F t (x l ,t)}-E T (x i ,t) + E s (x l ,t), 



where E$(xi, t) is the inverse Laplace transform of the error E${xi, s). 
Approximation errors Et and Ep 

The error Et that comes from the integral approximation using the trapezoidal rule, 
according to Crump [6], is 



E T = J2 e~ 2n ^ T n(x u 2nT + t). 



n=l 



Assume now that our function is bounded such as |n(xj,t)| < e CTi , for all Xj. Therefore the 
error can be bounded by 



E T < e- £ e-(^) = _ p < t < 2T. 

n=l 
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It follows that by choosing (3 sufficiently larger than a, we can make Et as small as desired. 
For practical purposes and in order to choose a convenient /3 we use the inequality which 
bounds the error 

£ T < e <rf-2T(/3-a)_ 

If we want to have the bound Et < 6t then by applying the logarithm in both sides of 
the previous inequality we have 

IT + 1 1 



Assuming a > we can write 



P > <r 



In i 



2T 



In our example we consider a = 0. In practice the trapezoidal error Et is controlled by 
the parameter f3 we choose. 

The second error, Ep 1 , comes from the approximation of the continued fraction given 
by (I26p . This error is controlled by imposing a tolerance TOL such as 

\v(z,M) -v (z,M-l)\ <TOL, 

in order to get the approximation rji(t) given by 



V. 



i(t) = ^Re{« (z, Mi)}, 



where Mi changes according to which Xi we are considering. 
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Figure 2: Error £ F and E G for P = 2, t = 1, < x < 12 and /3 = -ln(l(T 6 )/2T and 
different values of TOL. The global error is controlled by the parameter j3. 

In order to understand better how to control the trapezoidal error with the parameter 
/3 and how the tolerance TOL affects the error, we present a test example which is an 
analytically exactly solvable model. We assume P constant and Fickian diffusion 



dn 
~dt 



n 



p d_n +D &_ 

dx dx 2 



x e ]0,oo[,t > 0. 



The initial condition is n(x, 0) = 0, and the boundary conditions are 

n(0,t) = N , n(oo,i) = 0. 



(29) 



(30) 
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Figure 3: Error Ep and E G for P = 2, t = 1, < x < 12 and /? = - ln(lCT 10 )/2T 
and different values of TOL. The parameter /3 is chosen such that the global error is not 
affected. 



It will be noted that we are now considering a semi-infinite geometry. We note the dif- 
ference between this test case and our original unbound problem. We choose this test 
example for two reasons: Firstly, equation (I29p can be analytically exactly solved by first 
applying the time-Laplace transform and then through the inverse Laplace transform. Sec- 
ondly, this example is chosen to compare the convergence aspects of the Laplace inversion 
algorithm without spatial discretization. 

If we apply the Laplace transform to this problem we obtain 



n(x, s) = Na—e 
s 

The analytical solution is given by 



No 
2 



1 P/2D-x^/{P/2D)i+s 



(31) 



^erfc 


~x - Pt 




_2y/Di _ 



+ e p ^erfc 



x + Pt 



2VDt 

In Figures [5] and O for P = 2, t = 1 and < x < 12, we plot the following errors 



(32) 



Ep 



max 

Ki<N- 



\v(z,Mi)-v(z,Mi-l)\ 



and 



E G = \\n(xi,t) - r)i(t)\\ c 



(33) 



(34) 



where || • ||oo is the infinity norm. We choose the interval < x < 12 in order to avoid the 
influence of the right numerical boundary condition in the numerical computations, that 
in this case is n(12,t) = 0. 

The error Ep is related with the error Ep 1 since we control E^ by controlling Ep 
with the tolerance TOL. Figures [2] and [3] show how the parameter j3, given by (5 = 
- ln(lCT 6 )/2T in Figure [2] and /3 = - ln(10" 10 )/2T in Figure O affects the global conver- 
gence. Note that in Figure [2] the precision does not go further than 10 -6 . The global error 
of Figure [2] and Figure [3] is not affected by the spatial error E$ since we apply the Laplace 
inversion algorithm directly in (|3ip . 

The Laplace inversion algorithm approximates the value of the infinite series using a 
truncated continued fraction and this truncation is done by choosing an M, for each Zj, 
This Mi is chosen according to which value of the tolerance TOL we consider. We show 
in Figure H]the variations of Mi and it is clear the algorithm concentrates the high values 
of M in the region that presents steep gradients. 
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Figure 4: Number of iterations M for P = 2, t = 1 and TOL = 1/N 2 . Left figure: 
Approximate solution; Right figure: Number of iterations for each Xj. 



Spatial discretization error E$(xi, s) 

We now turn to the discretization error Es(xi,s), defined in (|'27|) (our main prob- 
lem), and prove that the method has truncation error of second order. Let us denote the 
differential operator L given by 

d 2 n d 
m = l^- X ^-d-x {P ^- 

We also denote by L n , the operator associated with the spatial discretization, given by 

L v n{xi,s) = — 2 Xsni(s) — , 

where rtj(s) denotes the exact solution at n(xi,s). The local truncation error is given by 

T e (xi, s) = L n n(xi, s) - Ln(xi, s). 

For a fixed s, we make a Taylor expansions of the functions n and P around the point Xj. 
We obtain, for a sufficiently smooth n, 

iTj i j_i(s)ni_i(s) + K iyi (s)ni(s) + Ki^ +1 (s)n i+1 (s) + (1 + s)n(xi, 0) 
d 2 rij . . . _ , . d 



(s) - A s ra;(s) - —(Pn)i(s) + (1 + s)n(x;, 0) 
ax 

1 ,„ \ n &i 1 ,d 2 ni l v d 3 ni 1 d 4 n; 

77 P n<(s) - 7;Pi -t-(s) ~ ~Pi-rr{s) ~ 77^ "TIT ( s ) + 



Ax 2 



6 1 w 2 1 dx y 1 2 1 dx 2 w 6 dx 3 w 12 dx 4 
+0(Ax 3 ), 

where P', P", P'" denotes the derivatives of P (unlike in the previous test example P is 
now not a constant). From this result we can conclude that, for n(., s) E C 4 (M), we have 

1 1 Xk 1 1 oo = max |T e (xi, s)\ < cAx 2 . 

2<i<JV 

By denoting Pj = Ps(xj, s), i = 1, . . . , JV — 1 we have 

L-n-Ei — T' e (xj, s), 
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that is, 

K(s)E{s)=T e (s). 

If ||ET _1 (s)|| 00 < C then \Ei\ < CHTeH^. Since the matrix K{s) is not an M-matrix 
|18| I19j. it is not easy to prove analytically the inverse of K{s) is bounded. This difficulty 
is related to the set of values of the parameter X s , given by 



A., 



+ a = P 2 + - u? + iu(2fi + 1) 



kir 



k = 1, 



where r defines the set of values in the Laplace domain, since for co 2 > j3 2 + /3 the complex 
A s has negative real part. However, it is easy to see numerically that for a fixed T, where 
T defines the stepsize of the trapezoidal rule used to approximate the integral (|23|) . as we 
refine the space step, the value ||-fC _1 (s)|| 00 does not change significantly. We also observe 
that ||if _1 ( s )ll 

oo is larger for values of \s\ close to zero, indicating that the convergence 
can be slower for these values, as can be observed in Figure [5j 



2 
|S| 



0.1 



0.2 
|s| 



0.3 



0.4 



Figure 5: Infinity norm of the matrix K 1 (s) for N = 1000. We have considered P(x) 
with a = 1. Left figure: T = 80; Right figure: T = 800. 

Additionally we observe that we have a similar phenomenon! to the so-called pollution 
effect p] observed for the Helmholtz equation and high wavenumbers where the discretiza- 
tion space step has to be sufficiently refined to avoid numerical dispersion. Also in this 
context it is observed that if we have a complex number as a coefficient in the equation, 
which is our case with X s , the imaginary part acts as an absorption parameter, which seems 
to allow us to control better the solution, decreasing the solution magnitude [9]. Following 
what is reported in literature [3l [16], a natural rule observed for an adjustment of the 
space step is to force some relation between T and the Ax. In our numerical experiments 
in order to retain accuracy we have considered 



* a 2 

-Ax < . 

T ~ 100 



(35) 



Numerical tests: Order of convergence of the numerical method 

In order to illustrate the behaviour of the numerical method, we consider two test 
problems. First, we consider the Fickian problem (I29p - (|30p . which exact solution is given 
by (|32p . In Table 1 we show the global error (|34p . for P = 2, t = 1 and different values of 
the space step. The value of T = 20 was considered in order to verify (I35p . These results 
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Ax 


Error 




Rate 


10/32 


0.4000 x 10" 


-2 




10/128 


0.2596 x 10" 


-3 


2.0 


10/512 


0.1625 x 10" 


-4 


2.0 


10/2048 


0.1016 x 10" 


-5 


2.0 



Table 1: Fickian case: Global error JMD for P = 2, t = 1, < x < 10, TOL = 1/iV 3 , 
/3 = -ln(10~ 16 )/2T, T = 20. 



confirm that the convergence is of second order. We can obtain similar results for other 
values of P. 

We now consider a non-Fickian problem given by the telegraph equation, 

f? + !? = D §- * £ ]0,2.M>0, (36) 

with initial conditions 

x Qtl I x 

n(x, 0) = sin(-), — (x, 0) = -- sin(-), (37) 

and boundary conditions 

n(0,t)=0, n(2vr,t)=0. (38) 
We can easily obtain the analytical solution given by 

n(x,t) = e~isin(|). (39) 

As for the Fickian case, we present in Table 2 the global error ()34j) for t = 1 and different 
space steps. We use the same value of T in order to have (f35|) . 



Ax 


Error 




Rate 


2?r/32 


0.5733 x 10" 


-4 




2vr/128 


0.4060 x 10" 


-5 


1.9 


2vr/512 


0.2397 x 10" 


-6 


2.0 


2vr/2048 


0.1627 x 10" 


-7 


1.9 



Table 2: Non-Fickian case: Global error (JM]) for t = 1, < x < 2vr, TOL = 1/N 3 , 
p = -ln(10- 16 )/2T, T = 20. 



It will be noted that our main problem is unbounded. But with at least one zero 
boundary condition for each, the two test examples ( Fickian and non-Fickian), although 
semi-bounded and bounded respectively, can be computationally viewed as similar to our 
main problem. We observe from Tables H] and H] for the test examples that we obtain 
second order convergence as predicted by the theoretical analysis for the main problem. 

5 Numerical results for n(x,t), j(x,t) and < x 2 (t) > 

To do the numerical experiments we consider the equations 

d 2 n dn d . d 2 n , An . 

W + m=-8-J Pn) + d^ (40) 



n 



and 



for 



P(x) 



dV 

dx 



dj dn 



and V(x; a) 



(41) 



1 



Jo(ia) 




40 



2(1 



-20 



Figure 6: Potential -P(x). Left figure: a = 1; Right figure: a = 16. 

For a = 1, the potential is smoother compared with a = 16, as it is shown in Figure 
[6l In Figure [6] we observe that P{x) for a = 1 changes between —1 and 1, whereas for 
a = 16 changes between —20 and 20 and the change is not smooth. Our method can deal 
very well with both cases. 

We consider the initial conditions, 



n (x, 0) 
j(x,Q) 



1 



-x 2 /L 2 



-x 2 /L 2 



£co) = o, 



1 

V- 



2x + P(x] 



(42) 
(43) 



and the boundary conditions are given by 



lim n(x,t) = and lim n(x,t) = 0. 

x— >oo X— >— oo 

Note that the stationary solution of the problem is given by 

n st (x) = N r exp (-V(x)) , 

where iV r is a normalization value. 

For a = 1 we show in Figures [7] and [8] the behaviour of the solution as we increase time 
from t = 1 until t = 500. The peak starts to split into two and then we have several waves 
forming that goes to the right and left. The domain where the function is not zero becomes 
larger as we travel in time. For that reason the computational domain increases consider- 
ably which requires more computational effort regarding the discretization in space. For 
an iterative method where we need to consider a discretization in time, it would require 
more computational effort for long times as we need to iterate in time whereas the Laplace 
transform has the advantage of not iterating in time and therefore it is the same if we 
compute the solution for short times or long times. 
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Figure 8: Particle density n(x,t) for a = 1. Left figure: Curve for instant of time t = 100 
;Right figure: Curve for instant of time t = 500. 



In Figure [9] we plot the flux for a = 1 , as it evolves from t = to t = 1 and in Figure 
[TOl as it evolves from t = 5 to t = 30. 

A quantity of physical interest in diffusion problems is the mean square displacement 
defined by 

/oo 
[x 2 n(x,t)]dx. 
-oo 

For the Fickian case, < x 2 (t) > is linear in t for all times in the absence of a potential. 
Now we would like to present calculations of < x 2 (t) > for the non-Fickian diffusion. At 
short times, and in the presence of a potential, the mean square displacement, < x 2 (t) >, 
shows a t 2 behaviour, see Figure 11. This is due to inertial effects which are captured by 
a non-Fickian diffusion equation. 

For a = 16 we show the evolution of the solution in the first instants of time. We 
see the solution presents very steep gradients and the method is able to give accurate 
solutions. First we observe how the wave split for t = 1 and t = 2 in Figure [121 

Next in Figure [13] we observe the behaviour for very large times. It is interesting to 
see how the Laplace method is able to give very quickly solutions for very large times. 
An iterative numerical method in time, it would take a large amount of time to run 
experiments for such long times such as t = 5000 or t = 10000 as we can see in Figure [Til 
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Figure 9: Density flux j(x,t) for a = 1. Left figure: Curve for instant of time t = 0; Right 
figure: Curve for instant of time t = 1. 





Figure 10: Density flux j(x,t) for a = 1. Left figure: Curve for instant of time t 
Right figure: Curve for instant of time t = 30. 



The flux for a = 16 is plotted in Figure [T5l 

6 Summary and outlook 

In this paper we have presented a numerical solution of a non-Fickian diffusion equation 
which is a partial differential equation of the hyperbolic type. This equation is of phys- 
ical interest in the context of Brownian motion in inertial as well as diffusive regimes. 
In our model the Brownian particle is subjected to a symmetric periodic potential of 
flexible shapes (generated with a single variable parameter) which can lead to harmonic, 
anharmonic or a confining potential for the particle. 

Instead of introducing discretization in both space and time variables we dealt with 
the time-derivatives through time Laplace transform and obtained an ordinary differential 
equation in space variable. This equation was then solved with a finite-difference scheme, 
leading to a discretised approximate solution for n(x, s); the solution is approximate due to 
discretization and is still formally exact in the Laplace domain s. The next step consisted 
of numerical Laplace inversion to obtain an approximation to the original solution n(x, t). 
We show that the full method is second order accurate; this finding receives additional 
support from two test examples considered in Section 4. One may be able to consider 
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Figure 11: Mean square displacement for a = 1; Left figure: Curves for t £ [0,2]; Right 
figure: Curves for i £ [5,30]. 
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Figure 12: Particle density n(x,t) for a = 16. Left figure: Curve for instant of time t = 1; 
Right figure: Curve for instant of time t = 2. 

further improvement. A major advantage of using the time Laplace transform is that we 
can compute the approximate solution for long times accurately and quickly. Any iterative 
numerical method would take too long to compute the solution for similar times even if we 
consider an unconditionally implicit numerical method which will allow large time steps. 
Additionally, our algorithm takes into consideration the smoothness of the solution; in 
other words the computational effort is higher in the regions where the solution has steep 
gradients. Another merit of the method is that it can be easily generalized to higher spatial 
dimensions. It would be of interest to consider an application of the method to numerically 
solve the Kramers equation which is a more involved partial differential equation than the 
non-Fickian diffusion equation. 
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